This is the analysis of all the salinity and nutrient experiments conducted on the lanai in St. John 616 from September 2021 to October 2022. These experiments incorporated four paired salinity and nutrient treatments with three temperatures. Each of the first four runs produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled in February, April, and October 2022. This output reflects all data totaling five treatments for Ulva and six treatments for Hypnea.
Packages loaded:
library(lme4)
library(lmerTest)
library(afex)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)
Load this file for normalized to quantum efficiency of photosynthesis per same as above
all_runs_photosyn_data <- read.csv("../data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")
Assign numerous variables as factors
all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)
all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)
all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))
all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)
Subset data to get only Ulva data for day 9 Ek– also remove the treatment 2.5 that was problematic
ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9 & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol"
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol"
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol"
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"
Ulva Ek________________________________________________________________________________________
hist(ulva$ek.1, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)
ulva %>% ggplot(aes(ek.1)) +
geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
Ulva Ek model
ulva_ek_model <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)
Make residual plots of the data for ulva ek
hist(resid(ulva_ek_model))
plot(resid(ulva_ek_model) ~ fitted(ulva_ek_model))
qqnorm(resid(ulva_ek_model))
qqline(resid(ulva_ek_model))
Check the performance of the model, get r2, make a table, plot effects, summary
performance::check_model(ulva_ek_model)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(ulva_ek_model)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## R2m R2c
## [1,] 0.6248358 0.7151076
summary(ulva_ek_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
## (1 | RLC.Order)
## Data: ulva
##
## REML criterion at convergence: 2145.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.75984 -0.58531 0.07584 0.56071 2.65765
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 76.53 8.748
## Run (Intercept) 43.24 6.576
## RLC.Order (Intercept) 18.86 4.343
## Residual 437.52 20.917
## Number of obs: 240, groups: Plant.ID, 95; Run, 7; RLC.Order, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 34.248 6.522 8.394 5.251 0.000661 ***
## Treatment1 39.125 7.213 6.771 5.425 0.001097 **
## Treatment2 46.727 7.213 6.771 6.479 0.000392 ***
## Treatment3 79.171 7.213 6.771 10.977 1.47e-05 ***
## Treatment4 84.260 7.213 6.771 11.682 9.80e-06 ***
## Temperature27 -10.542 4.842 32.140 -2.177 0.036904 *
## Temperature30 -10.619 4.435 63.820 -2.394 0.019600 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.688
## Treatment2 -0.688 0.825
## Treatment3 -0.688 0.825 0.825
## Treatment4 -0.688 0.825 0.825 0.825
## Temperatr27 -0.355 0.002 0.002 0.002 0.002
## Temperatr30 -0.344 0.000 0.000 0.000 0.000 0.475
tab_model(ulva_ek_model)
| ek.1 | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 34.25 | 21.40 – 47.10 | <0.001 |
| Treatment [1] | 39.12 | 24.91 – 53.34 | <0.001 |
| Treatment [2] | 46.73 | 32.52 – 60.94 | <0.001 |
| Treatment [3] | 79.17 | 64.96 – 93.38 | <0.001 |
| Treatment [4] | 84.26 | 70.05 – 98.47 | <0.001 |
| Temperature [27] | -10.54 | -20.08 – -1.00 | 0.030 |
| Temperature [30] | -10.62 | -19.36 – -1.88 | 0.017 |
| Random Effects | |||
| σ2 | 437.52 | ||
| τ00 Plant.ID | 76.53 | ||
| τ00 Run | 43.24 | ||
| τ00 RLC.Order | 18.86 | ||
| ICC | 0.24 | ||
| N Run | 7 | ||
| N Plant.ID | 95 | ||
| N RLC.Order | 6 | ||
| Observations | 240 | ||
| Marginal R2 / Conditional R2 | 0.625 / 0.715 | ||
plot(allEffects(ulva_ek_model))
Construct null model to perform likelihood ratio test REML must be FALSE
ulva_ek_treatment_null <- lmer(formula = ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
## boundary (singular) fit: see help('isSingular')
ulva_ek_model2 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_ek_treatment_null, ulva_ek_model2)
## Data: ulva
## Models:
## ulva_ek_treatment_null: ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_ek_model2: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## ulva_ek_treatment_null 7 2332.8 2357.1 -1159.4 2318.8
## ulva_ek_model2 11 2200.1 2238.3 -1089.0 2178.1 140.69 4 < 2.2e-16
##
## ulva_ek_treatment_null
## ulva_ek_model2 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_ek_temperature_null <- lmer(formula = ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_ek_model3 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_ek_temperature_null, ulva_ek_model3)
## Data: ulva
## Models:
## ulva_ek_temperature_null: ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_ek_model3: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## ulva_ek_temperature_null 9 2202.8 2234.2 -1092.4 2184.8
## ulva_ek_model3 11 2200.1 2238.3 -1089.0 2178.1 6.7989 2
## Pr(>Chisq)
## ulva_ek_temperature_null
## ulva_ek_model3 0.03339 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Plot raw data
ulva %>% ggplot(aes(treatment_graph, ek.1)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="salinity/nitrate", y= "Day 9 Ek (μmols photons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-1, 250) + stat_mean() +
scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
Make linear regression with growth add growth rate from other dataset to this one and subset by species
growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$lunar.phase <- as.factor(growth_rate$lunar.phase)
Make a new column for weight change (difference final from initial) add growth column to dataset
growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
ulva$lunar.phase <- (gr_ulva$lunar.phase)
Plot linear regression between ek and growth
ulva_growth_ek_graph <- ggplot(ulva, aes(x=ek.1, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Ulva lactuca Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)",
y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_ek_graph
## `geom_smooth()` using formula = 'y ~ x'
Summarize the means for Ek
ulva %>% group_by(Treatment) %>% summarise_at(vars(ek.1), list(mean = mean))
## # A tibble: 5 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 27.2
## 2 1 65.2
## 3 2 72.8
## 4 3 105.
## 5 4 110.
Hypnea Ek____________________________________________________________________ Use Day 9 for final analysis hm1-1 on 10/12/22 and hm1-2 on 4/29/22 causing issues of influential observations hm1-2 for both rETRmax and Ek – leaving them in dataset because no good reason to believe not good data There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21 below to match growth data
hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol"
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol"
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol"
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"
Make a histogram for hypnea
hist(hypnea$ek.1, main = paste("Hypnea musciformis -- Ek"), col = "maroon2", labels = TRUE)
hypnea %>% ggplot(aes(ek.1)) +
geom_histogram(binwidth=5, fill = "#7D0033", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Run model for Hypnea and Ek
hyp_ek_model <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)
Plot residuals
hist(resid(ulva_ek_model))
plot(resid(hyp_ek_model) ~ fitted(hyp_ek_model))
qqnorm(resid(hyp_ek_model))
qqline(resid(hyp_ek_model))
Check the performance of the model, get r2, make a table and plot for model
performance::check_model(hyp_ek_model)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(hyp_ek_model)
## R2m R2c
## [1,] 0.3763202 0.5670999
summary(hyp_ek_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
## Data: hypnea
##
## REML criterion at convergence: 2625.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0609 -0.6195 -0.0476 0.4318 3.7844
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 112.8 10.62
## Run (Intercept) 127.8 11.31
## Residual 546.1 23.37
## Number of obs: 286, groups: Plant.ID, 143; Run, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 76.847 9.143 6.026 8.405 0.000151 ***
## Treatment1 18.182 10.840 6.057 1.677 0.144009
## Treatment2 17.509 10.840 6.057 1.615 0.156904
## Treatment2.5 42.029 14.806 4.602 2.839 0.039907 *
## Treatment3 17.303 10.840 6.057 1.596 0.161073
## Treatment4 24.523 10.864 6.113 2.257 0.063977 .
## Temperature27 -35.829 4.237 102.763 -8.456 1.98e-13 ***
## Temperature30 -39.418 4.224 106.633 -9.331 1.71e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.784
## Treatment2 -0.784 0.903
## Treatmnt2.5 -0.574 0.484 0.484
## Treatment3 -0.784 0.903 0.903 0.484
## Treatment4 -0.783 0.901 0.901 0.483 0.901
## Temperatr27 -0.231 0.002 0.002 0.000 0.002 0.002
## Temperatr30 -0.231 0.000 0.000 0.000 0.000 0.005 0.496
tab_model(hyp_ek_model)
| ek.1 | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 76.85 | 58.85 – 94.85 | <0.001 |
| Treatment [1] | 18.18 | -3.16 – 39.52 | 0.095 |
| Treatment [2] | 17.51 | -3.83 – 38.85 | 0.107 |
| Treatment [2.5] | 42.03 | 12.88 – 71.18 | 0.005 |
| Treatment [3] | 17.30 | -4.04 – 38.64 | 0.112 |
| Treatment [4] | 24.52 | 3.14 – 45.91 | 0.025 |
| Temperature [27] | -35.83 | -44.17 – -27.49 | <0.001 |
| Temperature [30] | -39.42 | -47.73 – -31.10 | <0.001 |
| Random Effects | |||
| σ2 | 546.12 | ||
| τ00 Plant.ID | 112.83 | ||
| τ00 Run | 127.85 | ||
| ICC | 0.31 | ||
| N Run | 8 | ||
| N Plant.ID | 143 | ||
| Observations | 286 | ||
| Marginal R2 / Conditional R2 | 0.376 / 0.567 | ||
plot(allEffects(hyp_ek_model))
Construct null model to perform likelihood ratio test REML must be FALSE
hyp_ek_treatment_null <- lmer(formula = ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_ek_model2 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_ek_treatment_null, hyp_ek_model2)
## Data: hypnea
## Models:
## hyp_ek_treatment_null: ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_ek_model2: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## hyp_ek_treatment_null 7 2690.8 2716.3 -1338.4 2676.8
## hyp_ek_model2 12 2690.0 2733.8 -1333.0 2666.0 10.789 5 0.05574
##
## hyp_ek_treatment_null
## hyp_ek_model2 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hyp_ek_temp_null <- lmer(formula = ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_ek_model3 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_ek_temp_null, hyp_ek_model3)
## Data: hypnea
## Models:
## hyp_ek_temp_null: ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_ek_model3: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## hyp_ek_temp_null 10 2740.4 2777.0 -1360.2 2720.4
## hyp_ek_model3 12 2690.0 2733.8 -1333.0 2666.0 54.454 2 1.498e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Use ggplot for the data in a boxplot
hypnea %>% ggplot(aes(treatment_graph, ek.1)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="salinity/nitrate", y= "Ek (μmols photons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-1, 250) + stat_mean() +
scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
For Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/9/21 because it was white and also looked dead
gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)
Plot a regression between the photosynthetic independent variables of interest and growth rate
hypnea_growth_ek_graph <- ggplot(hypnea, aes(x=ek.1, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Hypnea musciformis Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)",
y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor(label.y = 150)
hypnea_growth_ek_graph
## `geom_smooth()` using formula = 'y ~ x'
Summarize the means
hypnea %>% group_by(Treatment) %>% summarise_at(vars(ek.1), list(mean = mean))
## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 51.8
## 2 1 71.4
## 3 2 70.7
## 4 2.5 93.8
## 5 3 70.5
## 6 4 78.3
hypnea %>% group_by(Run) %>% summarise_at(vars(ek.1), list(mean = mean))
## # A tibble: 8 × 2
## Run mean
## <fct> <dbl>
## 1 1 90.1
## 2 2 79.1
## 3 3 60.8
## 4 3.5 64.5
## 5 4 58.7
## 6 7 93.8
## 7 8 52.2
## 8 9 51.3